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ABSTRACT 

o . 

f^ , We investigate the distribution of mass M and orbital period P of extra-solar planets, 

^^ \ taking account of selection effects due to the limited velocity precision and duration 

Q ' of existing surveys. We fit the data on 63 planets to a power-law distribution of the 

O ; form dn = CM-°'P-f^{dM/M){dP/P), and find a = 0.12 ± 0.10, /3 = -0.26 ± 0.06 

p\J ' for M ^ lOMj, where Mj is the Jupiter mass. The correlation coefficient between 
these two exponents is —0.32, indicating that uncertainties in the two distributions 

p\J . are coupled. We estimate that 3% of solar-type stars have companions in the range 

> ; IMj < M < lOMj, 2 d < P < 10 yr. 

OO 

f^ I 1. Introduction 



o 



S 



As of May 2001, radial- velocity surveys have discovered over sixty planets orbiting nearby stars. 
This sample should be large enough to provide reliable estimates of the distribution of planetary 
mass and orbital elements, at least in the range to which the radial-velocity surveys are sensitive. 



At least two important selection effects must be included in any statistical analysis of this 
C^ ■ kind: (i) each survey has a detection limit K^, such that the orbits of companions that induce 

refiex motions in their host star of amplitude < Kd cannot be reliably characterized; (ii) orbits of 
^ , companions with periods much longer than the duration of the survey cannot be reliably character- 

ized. In any survey limited by its velocity precision, uncertainties in the distribution of planetary 
masses M are coupled to uncertainties in the distribution of orbital periods P, because the velocity 
amplitude induced by a companion is oc MP~^'^ (eq. ^). Thus it is necessary to determine both 
distributions simultaneously. 

The aim of this paper is to describe a maximum-likelihood method of estimating these distribu- 
tions using data from a variety of surveys, while accounting for survey-dependent selection effects. 
We have chosen to fit the data to simple power-law models of the distribution of masses and peri- 
ods; such distributions are simple to interpret and common in nature, and it is straightforward to 
generalize our approach to non-parametric models as the data improve. 

We shall work with data from eight radial- velocity surveys of nearby stars, which have detected 
between 2 and 23 extra-solar planets each (see Table |l|). Together these surveys have detected 90 
planets as of June 1 2001, although several planets appear in more than one survey so we have only 
63 distinct planets. 



We compare our method and results with other determinations of the mass distribution of 
extra-solar planets in §0. 



2. Maximum Likelihood Method 

We focus initially on the simple case of a single survey that examines A^* stars for radial velocity 
variations due to orbiting companions. The velocity amplitude K due to a companion of mass M 
with orbital period P is 

where e is the orbital eccentricity, M^ is the stellar mass, and i is the inclination between the orbit 
plane and the sky plane. For simplicity we assume that e <C 1 so that the factor (1 — e^)~^'^ 
is unity. We do not attempt to account more accurately for the eccentricity dependence because 
the detectability limit for eccentric orbits depends on both the amplitude and shape of the radial- 
velocity curve. At the upper quartile of the eccentricities in our sample, e = 0.46, the error in 
K caused by setting e = in equation (|l]) is only 11%. We also assume that M <C M*, so that 
equation (|l|) simplifies to 



K = ^ 

M* 



27rGM* 



1/3 

p , (2) 

where m = M sin i is the minimum companion mass, corresponding to an orbit viewed edge-on. 

Throughout this paper we shall assume that all of the stars in the survey have mass equal 
to the Sun's, M-i, = Mq. This is not a bad approximation since most radial- velocity surveys for 
low-mass companions have focused on solar- type stars. In fact, our method is easily generalized to 
the case where the survey stars have different masses, but to do so we need to know the masses 
of all the stars in the survey (not just the ones that have detected companions) — and this extra 
complication did not seem worthwhile in this preliminary analysis.Q 

We restrict our attention to companions with minimum mass m < mmax = lOMj, where Mj is 
the Jupiter mass; this cutoff hopefully minimizes the contamination of our sample by brown dwarfs 
and is below the deuterium-burning threshold which sometimes is taken to define the boundary 
between planets and stars. We also restrict our attention to orbital periods P > Pmm = 2 days, 
corresponding to a semimajor axis of 6.7Rq = 0.031 AU; this limit is small enough to include all 
the known planets. 

We assume that the probability that a single star has a companion with mass and orbital 



^One of the largest relative errors caused by setting e = and M, = AIq is for e Eri (Af* — O.SMq, e = 0.61), 
where equation (Q) with Af, — Mq yields a value for K that is 46% too small. 



period in the range [M, M + dM], [P, P + dP] is given by a power law, 

/ M \ -" / P \ -^ dM dP 



where C, a and /3 are constants to be determined, and Mq and Pq are a fiducial mass and period, 
which we choose to be Mq = 1.5Mj and Pq = 80 d (the reasons for this choice are outlined in 
the following subsection). If the distribution of companion orbits is isotropic, the distribution of 
minimum mass m = Msini and period is given by 

f mX-'^ f P\-l^ dmdP 



where 



,r(i + i 



a 



,2 



-C, (5) 



r(2 + Q) 

r(-) is the gamma function, and q > — 2. 

Initially we assume that the survey detects a companion if and only if (i) the velocity amplitude 
K exceeds a survey-dependent detectability limit Kd; and (ii) its orbital period is shorter than a 
survey-dependent upper limit Pmax- We expect that Pmax will be proportional to the duration of 
the survey, since typically at least two orbits are required for a reliable detection. More realistic 



smooth cutoffs to the detection efficiency are discussed in § |2.2| . Although the detection limits Kd 
and Pmax can be estimated from descriptions of the survey, we adopt the more objective approach 
of determining them directly from the maximum-likelihood analysis. 

Let Xi = ln(mj/Mo) and yi = ln(Pj/Po), i = 1, . . . , A^, where m,j and Pj are the minimum 
mass and period of the companions detected in the survey. Then the velocity amplitude Ki (eq. |2|) 
exceeds the detection limit Kq if 

Similarly, the orbital period is less than the maximum detectable period if 

y.<^ = ln(^). (7) 

The other constraints are 

Xi < Xmax = ln(mmax/Mo) = In 10 + ln(Mj/Mo), 

Vi > ymm = hi(Pmin/Po) = -5.207 + ln(lyr/Po). (8) 

The constants Xmax and ymin are fixed, while the variables u and v are to be determined by the 
maximum-likelihood analysis. 

From equation (Q) the expected number of companions in the interval dx dy in a survey of N* 
stars is 

n{x,y)dxdy = N*p{x,y)dxdy where p{x,y) = ce~"^~^^. (9) 



The likelihood function L is the product of (i) the probability of detecting N companions with 
minimum masses Xj and periods y^; and (ii) the probability of observing none elsewhere in the 
domain D of (x, y) space in which companions are detectable. Thus 



TV . 

L = TTn(xj,2/i)exp -/ dxdyn{x,y) 
t=i L Jd 



if aU {xi,yi) £ D, 



and zero otherwise. The domain D is v + ^y < x < Xmax, ymin < y < u, with 

u{u, v) = min [u, 3(a;max - v)] . 
We now substitute equation (pi) into equation ^^) and take the log of the result, 



N 



N 



In L = TV ln(cA^*) -a^Xi-/?^yi- cN*f{a, (3, u, v). 



i=l 



i=l 



Here 



where 



f{a,(3,u,v) 
g{a,P,x,y) = \ 



dy / dxg{a,P,x,y), 

— oo Jv+y/3 



^ ax Py jf ^ ^ ^^^^ ^^^ y y y^.^^ 

otherwise. 



The integral yields 



f{a,f3,u,v) 



3e- 



a{a + 3/5) 



e 3 



|(a+3/3)i^min _ „-|(a+3/3)u(M,t)) 
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(10) 

(11) 

(12) 

(13) 
(14) 

(15) 



if f < Xinax — -^ymin and u > i/min, and zero otherwise. 



The best estimates for the fitted variables c, a, (3, u and v correspond to the global maximum 
of InL. First, the constant c can be evaluated from 



oc c 

Substituting this result into equation (|l2|) yields 

lnL = iV ln( — 1 -1 






N N 

a^Xi- P^yi 

i=l i=l 



(16) 



(17) 



To determine the best estimates for u and v we note that InL depends on these parameters only 



through /(a, /?, u, v), and that InL is maximized when / is minimized. According to equation (13), 
/ depends on u and v only through the limits of integration, that is, only through the shape of the 
domain D. Since the integrand is non- negative, we minimize / by making D as small as possible, 
so long as it still contains all the data points {xi,yi). This can be done by setting v equal to the 
smallest value of Xi — |yj in the sample, and u equal to the largest value of yi in the sample. 



The best estimates for the remaining parameters are given by 



d\nL 
da 

dlnL 

dp 
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' f dp ^^' 



(18) 
(19) 



i=l 



which can easily be solved numerically. Table y lists the best estimates of the parameters for each 
survey. The value of the normalization parameters C is based on the fiducial mass and period 
Mq = 1.5Mj and Pq = 80 d, which are chosen for reasons outlined in the following subsection. 

Figure || shows the correlation between the duration of each survey and the fitted value of Pmax 
for that survey, as well as the stated velocity precision Ks for each survey and the fitted detection 
limit Ko for that survey (values taken from Tables |l| and ^). The detection limit is generally a 
factor of three or so higher than the stated precision, presumably because determining a reliable 
orbit is more difficult than simply detecting the presence of a companion. For most of the surveys 
there is an good correlation between duration and i-*max, with the slope of the correlation indicating 
that approximately two orbital periods of data are needed for a reliable detection. 



2.1. Generalization to multiple surveys 

It is straightforward to expand the analysis of the previous Section to multiple surveys, which we 
label by J = 1, . . . , J. The three parameters describing the companion distribution, a, /5 and c, are 
now derived from the entire sample of known companions from all surveys, while the parameters Uj 
and Vj that describe the period and radial-velocity thresholds are different for each survey. Thus 
the expected number of companions to be discovered in the interval dx dy in survey j is 



(20) 



nj(x, y)dx dy = N*p{x, y)dx dy, 



where p{x, y) is defined in equation (P) and N'- is the number of stars in survey j. Equation (OO 
becomes 



Af, 



Lj = ]^nj(xjj,2/ij)exp 



i=l 



D, 



nj{x,y)dxdy 



and the likelihood is given by 



In L = 2_. In Lj . 



(21) 



(22) 



Again, the integration constant can be eliminated from equation (22) 
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(23) 



and the likelihood becomes 



J 



J J 

lnAr;+ln^iV,-ln^iV;/,-l 
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As before, ijj is set equal to the smallest value of Xij — \yij in survey j, and Uj is set equal to the 
largest value of yi,j in survey j. 

The surveys listed in Table |l| have discovered 90 companions, although several appear in more 
than one survey so there are only 63 distinct companions. Companions discovered in multiple 
surveys are counted in each survey where they appear; this approach leads us to underestimate the 
statistical uncertainties in our parameters somewhat (probably by about a factor of (90/63)^'^ = 
1.2), but avoids the systematic bias that would be created by counting the companions only once 
and discarding them from the other surveys. A conservative alternative approach is to use only the 
Coralie survey for parameter estimation (top lines of Tables [l| and Q) . 

The values of the normalization parameters c and C quoted in this paper are based on the 
fiducial mass and period Mq = 1.5Mj and Pq = 80 d. These values are chosen to minimize the 
uncertainty in Inc. If the uncertainties are small, this requirement is equivalent to choosing Mq 
and Pq so that the covariances between c and the exponents a and f3 vanish. 

The likelihood analysis yields 

a = 0.12 ±0.10, (25) 

/3 = -0.26 ±0.06, (26) 

c = 1.7l|J:J? X 10"^ (27) 

C = 1.8t[J:J| X 10"^ (28) 

where the confidence limits correspond to InL = (InL)max— 0.5. With these estimators at hand, it is 
straightforward to plot the likelihood as a function of the exponents a and /3 (Figs. ^ and ^). Figure 
^ indicates that there is a significant covariance between the exponents a and j3 that characterize the 
mass and period distributions (correlation coefficient r = —0.32). This correlation arises because 
of the selection effects on velocity amplitude, and demonstrates that both distributions should be 
fitted simultaneously. 

Table y shows the difference between the number of predicted companions and the number of 
actual detections; these are generally in good agreement except for the McDonald survey, which is 



discussed further in S3.2 



This table also lists the number of predicted brown dwarfs (lOAfj < M < 80Mj) in each 
sample, assuming that the planetary mass function extends to 80Mj, and the number of brown 
dwarfs actually discovered. As many authors have pointed out, the small number of brown dwarf 
discoveries strongly suggests that the mass function we have derived cannot be extrapolated to 
brown dwarf masses. 



2.2. Generalization to a smooth cutoff 

A sharp cutoff in the detectabihty of planets at radial velocity K£) and period Pmax is not very 
realistic. A better approximation is to replace the sharp cutoffs in equations ( p!2| ) and ( [I^ ) with 
smooth functions. We can do this by replacing f{a,l3,u,v) with 

fs{a,P,u,v) = / dyhu{u-y) dxK{x - \y - v)e~°"'~'''y , 

•^J/min J -CO 

OO POO 

dyhuiu-y) dxK{x - \y - v)g{x,y), (29) 

'OO J — OO 

where g{x,y) is defined by equation (|j). 

The functions /lu(-) and h^{-) are measures of the detection efficiency of the survey as a function 
of orbital period and velocity amplitude. The function /iu(s) approaches as s ^ — oo and 1 as 
s — > oo; we shall assume that hu{s) — 2 is an odd function of s so that /iu(0) = 2, with similar 



assumptions for h^. In the limit where hu and h^ are step functions we recover equation (13). Thus 
V and u are still defined by equations (P) and (j^, except that Kd and Pmax are interpreted as the 
velocity amplitude and orbital period at which the detection efficiency falls to 50%. 

Let 

hu{s) = r bu{s')ds', Kit) = I K{t')dt'- (30) 

J —00 J— 00 

then equation (p9|) can be rewritten as 

/oo /"GO pU — s' poo 

ds'buis') / dt'h^{t') I dy dxg{x,y) 

-00 J —00 J —00 Jt'+y/3+v 

00 /-oo 

dsbu{s) / dtbv{t)f{a,(3,u — s,v + t), (31) 



where the second line follows from equation (|T3|), and we have dropped the primes on the dummy 
variables. These integrals are easy to evaluate numerically. 

We shall choose 

with a similar choice for hy{t). We call 6u and 6^ the threshold widths. Figure shows the effect of 
non-zero threshold widths on the slope estimators a and /5. In the arbitrary but plausible case where 
the detection efficiencies h drop from ^ to ^ over a factor of two in period or velocity amplitude, 
we have 5„, 6y = 0.51. In this case the best-fit value of a is shifted downward by 0.04 and the best 
fit for j3 is shifted upwards by about 0.020. These changes are significant but relatively modest; 
since the appropriate values of the threshold widths are difficult to estimate, we shall not attempt 
to correct for this effect. 



3. Discussion 



3.1. Comparison to other estimates of the mass distribution 

We have found that the distribution of companion masses below lOMj is approximately flat in 
logM, or slightly rising towards small masses, i.e., the exponent a is small and positive (a = 
0.12 =b 0.1, eq. p5|). The distribution of companion masses has already been examined by a number 



of authors, most of whom have reached similar conclusions (Mazeh et al. 1995, Marcy & Butler 
1998| , pVIazeh 199g| , ^tepinski fc Black 2000| , ptepinski fc Black 2001] [Jorissen et al. 200 1| , [Zuckei 



& Mazeh 2001). Our approach offers several advantages over the variety of methods used in these 
papers: (i) we correct for selection effects in period and velocity amplitude, (ii) we account for 
the coupling between the orbital period distribution and mass distribution, (iii) we estimate the 
sensitivity in radial velocity or maximum period (our parameters Kd and -Pmax) self-consistently 
from the data, and (iv) we determine the normalization of the mass distribution, not just its shape. 



3.2. Extrapolations 

It is interesting to investigate the implications of extrapolating the mass and period distributions 
that we have derived. If we assume that our maximum-likelihood distribution (eqs. 2^2^) applies 
in the mass range lOMj < m < 80Mj that is usually associated with brown dwarfs, we predict that 
the Coralie and Keck surveys should have discovered, respectively, 13 and 8 companions in this 
range (Table p[); in fact these two surveys found only one companion each. This result confirms 



the finding of several authors (Basri & Marcy 1997, Mayor et al. 1998, Mazeh et al. 1998, Mazeh 
1999, Jorissen et al. 2001) that there is a cutoff in the power-law distribution of companion masses 
at m ^ lOMj, and a "brown-dwarf desert" between ~ lOMj and ~ lOOMj in which few companions 
exist at semi-major axes less than a few AU. 

The average number of planets per star with masses between Mi and M2 and periods between 
Pi and P2 is given by equation (pi): 



N 



C_ 

a/3 



mJ 



Mo 
M2 



Pi. 



P2 



(33) 



Thus, for example, in our best-fit model (eqs. |25-28), the expected number of planets per star with 
periods between 2 days and 10 yr and masses between M and lOMj is 



N = 0.129 



(Mo/M) 



0.12 



0.755 



(34) 



For M = Mj, A^ = 0.033; thus about three percent of solar-type stars have a planet of Jupiter 
mass or larger in this period range. If we make the large extrapolation to Earth-mass planets 
(M = 0.003Mj) we find A^ = 0.172; in this case, more than 15% of stars would have an Earth-mass 
or larger companion. 



If we extrapolate to larger orbital periods, we find that the number of companions in a given 
mass range with periods between 2 days and 5 yr would be about 0.28 times the number with 
periods between 5 yr and 1000 yr; in this case a significant fraction of Jupiter-mass planets would 
have orbital periods short enough to be detected in existing radial- velocity surveys. 



3.3. Comparison with the solar nebula 

The mass and period distribution (^) can be used to determine the total surface density in planets 
less massive than Mmax; assuming that the central star has mass IMq: 

v^ ^ 3C Mq /M^axV"" ( Po Y (^ AU\2+3/3/2 

""^"^ = 4.(1 - a) (1 KW [-M^] [lYr] [-ir) ■ ^''^ 

For our best-fit model, 

S(a) = 50gcm (^^ (^--j ; (36) 

this can be compared to the gas and dust densities required in the minimum solar nebula (e.g. 

1 AU\i-^ 



Hayashi 19811) 



Sgas(a) = 1.7 X 10=^ g cm-2 



a 



1 AU\ ^'^ 
Sdust(a) = 7.1gcm-2( . (37) 



a 



The agreement of the exponents in equations (^) and (^) is striking and perhaps surprising, given 
that many theorists believe that the extrasolar giant planets must have formed at much larger radii 
and migrated to their present locations, while the planets in our solar system have suffered little 
or no migration. 



3.4. Summary 

Figure ^ shows the minimum-mass and period distributions of all of the substellar companions 
found in the surveys in Table |l|, along with the simple power-law models that we have used to 
fit these distributions. The effects of the selection effects at small mass and large period, and the 
evidence for a cutoff above ~ lOMj, are evident in the bottom panels. 

We have described a simple maximum-likelihood method that determines the mass and period 
distributions of extrasolar planets discovered in multiple surveys. Our method determines and 
accounts for selection effects on velocity amplitude and period from the data themselves, without 
relying on the nominal survey parameters. Our best-fit model is defined by equations (^) and 
equations (p5|)-(p^). 
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Survey Ks duration number of number of dis- references 

(m/s) (yr) stars observed covered planets 



Coralie 


4 
2-5 
3-10 


2.5 
5 
13 


1000 
530 
300 


23 
22 

17 


lUdry et al. 2001 


Keck 


Vogt et al. 200C 


Lick 


pumming et al. 1999 



Udry et al. 2001 



Nisenson et al. 199£ 



Tinney et al. 2001 



Endl et al. 2000 



Cochran et al. 2000 



Table 1: Characteristics of eight radial- velocity surveys of extrasolar planets: stated velocity pre- 
cision of the survey (Ks), duration of the program, number of stars observed to date, number of 
planetary companions discovered (Msini < lOMj), and references for each survey. 



Survey 


a 


/3 


CxlO^ 


-» max 

(yr) 


(m s" 


) 


Coralie 


0.16 ±0.19 


-0.23 ±0.13 


11-5111 


2.1 + 0.2 


11.1- 


0.7 


Keck 


0.06 ±0.19 


-0.15 ±0.12 


is.ilt? 


3.0 ±0.4 


10.3- 


0.9 


Lick 


-0.03 ± 0.23 


-0.09 ±0.12 


22.2l|9 


6.9 ±1.7 


12.9- 


1.5 


Elodie 


o.5il°:l? 


-n 44+0-17 


32.51^1,^ 


6.3 ±1.1 


36.6- 


1.9 


AFOE 


0.49t'oil 


-0.16l°:i 


re q+26.7 


3.6 ±3.3 


32.7- 


3.8 


AAT 


^■'^^-o'.m 


-i-^nti 


SI 4+19-8 


2.0 ±0.6 


38.6- 


2.8 


ESO 


0-98tlf, 


-, Q1+0.65 
-■^ -"J -'^-0.92 


A 0+15.8 
^•^-3.9 


6.9 ±2.0 


12.9- 


3.5 


McDonald 


2.09l\il 


-2A6l\f, 


n 1+2.5 
U-J--0.09 


6.9 ±1.0 


12.9- 


1.6 



Table 2: For each survey, best estimates for the exponents a and f3 and the normalizing constant 
C in equation (^), the maximum detectable period Pmax and the velocity precision Kd in m s~^. 
The fiducial mass and period are Mq = 1.5Mj, Pq = 80 d. The errors on Pmax and Kd are one- 
sided, since the maximum likelihood is achieved when these parameters equal the largest period 
and smallest velocity amplitude found in the survey. 
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Survey number of pre- number of dis- number of predic- number of discovered 

dieted planets covered planets ted brown dwarfs brown dwarfs 



Coralie 


351^ 


Keck 


2ltl 


Lick 


13 ±2 


Elodie 


9l? 


AFOE 


3±0 


AAT 


4±1 


ESO 


2±0 


McDonald 


o + l 
"J n 



23 

22 

17 

13 

7 

4 

2 

2 



13 



+7 



8i^ 

2±1 
3±1 
1±0 
1±1 



Table 3: Comparison of predicted and observed numbers of companions. 
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Fig. 1. — The left panel shows the correlation between the duration of each survey (from Table 
H) and the longest period in the sample. Note that Lick, ESO and McDonald programs have all 
detected the same long-period planet: e Eridani, which has an inferred orbital period of P = 2502.1 
days. All of the points lie below the dashed line {x = y), indicating that a reliable detection requires 
following the star for more than one orbital period. The right panel shows the correlation between 
the stated velocity precision of each survey (from Table |l|) and the smallest velocity amplitude of 
any of their detected planets. Almost all of the points lie well above the dashed line, indicating 
that determining a reliable orbit requires a velocity amplitude that is significantly larger than the 
stated velocity precision. The exception is the derived limit for McDonald; in this case K^ is set 
by their detection of a planet in the e Eri system, for which our approximations of a circular orbit 
and solar-mass star yield an estimate for the velocity amplitude that is 46% too low (see footnote 
1). 
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Fig. 2. — Contours of constant likelihood for the combination of all eight surveys. The maximum 
of L is located at the filled circle, q = 0.12 it 0.10, (3 = —0.26 it 0.06. The contour levels represent 
"n-cr" confidence regions, in which the likelihood function is smaller than its maximum value by 
exp(— n^/2). 
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Fig. 3. — The estimators of the exponents of the mass and period distribution, a and /?, for the 
combined eight surveys (fihed circles) compared with their values for each individual survey. The 
estimates are approximately consistent given the uncertainties. 
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Fig. 4. — The difference of the modified best parameters am and (3m (which depend on the threshold 



widths 5u and 6^) and the nominal values of a and (3 as quoted in eqs. (Esl) and (26). 
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Fig. 5. — The minimum mass (left panels) and period (right panels) distributions, on linear (top) 
and logarithmic (bottom) scales. The gray histograms represent all objects, including those with 
Msini > lOMj, found in the radial- velocity surveys in Table || (67 objects). The unshaded his- 
tograms include duplicate discoveries of the same object in different surveys, as discussed in § |2.1| 
(94 objects). The curves show the predictions of our best-fit model, given by equations (|25|)-Pq) 
and (^). As discussed in the paper, the histograms are subject to selection effects at small mass 
and large period, and there is evidence for a cutoff to the mass distribution above ~ lOMj. 



